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A new numerical algorithm for interacting fermion systems to treat the grand-canonical en- 
semble is proposed and examined by extending the path-integral renormalization group method. 
To treat the grand-canonical ensemble, the particle-hole transformation is applied to the Hamil- 
tonian and basis states. In the interaction-term projection, the Stratonovich-Hubbard trans- 
formation which hybridizes up and down spin electrons is introduced. By using this method, 
the phase diagram of the two-dimensional Hubbard model with next-nearest-neighbor transfer 
is accurately determined by treating the filling-control (FC) and bandwidth-control (BC) Mott 
transitions on the same ground. A V-shaped Mott insulating phase is obtained in the plane of 
the chemical potential and the Coulomb interaction, where the transitions at the corner (BC) 
and the edges (FC) show contrasted characters with large critical fluctuations near the edges 
coexisting with the first-order transition at the corner. This contrasted behavior is shown to 
be consistent with the V-shape structure of the phase boundary because of a general relation, 
in which the slope of the metal-insulator transition line in the phase diagram is expressed by 
thermodynamic quantities. The V-shaped opening of the Mott gap is favorably compared with 
the experimental results of the transition metal oxides. 

KEYWORDS: path integral renormalization group method (PIRG), grand-canonical ensemble, GPIRG, geometrical 
frustration, Mott transition, two-dimensional Hubbard model 



§1. Introduction 

The nature of the system in which quantum fluctua- 
tions and electron correlations play essential roles is one 
of the main subjects in condensed matter physics. When 
the kinetic energy and the Coulomb repulsion compete 
severely, the ground state of many-body electron sys- 
tems can be highly nontrivial. Metal-insulator transi- 
tions driven by the electron correlation postulated by 
Mott 1 - 1 provides a typical example of such a nontrivial 
behavior. 2 ) 

In the transition to the Mott insulator, it is known that 
there exist two different basic routes to control the com- 
petition of the interaction and kinetic energies. 2 ' One is 
the control by bandwidth (relative to the Coulomb re- 
pulsion) and the other is filling. Controls by these two 
parameters can be found in a lot of examples in real ma- 
terials including transition metal compounds, 2 - 1 organic 
materials 3 ' and 3 He systems. 4,5 ' In spite of plenty of 
the experimental results, phase diagrams of the Mott 
insulator and metals have not been fully elucidated in 
microscopic theoretical descriptions. 

The filling-control Mott transition (FCMT) was stud- 
ied at zero temperature by the quantum Monte Carlo 
(QMC) method 6 ' in the Hubbard model on a square 
lattice. 7 ' 8 ' The transition shows a continuous charac- 
ter with a singular divergence of the compressibility and 
critical divergence of the antiferromagnetic correlation 
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length. The bandwidth-control Mott transition (BCMT) 
was studied also at zero temperature by the path-integral 
renormalization group (PIRG) method. 3 ' 9 ' In contrast to 
the FCMT, the BCMT shows a first-order transition, al- 
though a naive expectation is that the continuous FCMT 
anticipates also the continuous BCMT. This contrast is 
essentially consistent with the trend of the experimental 
observations cited above. Mott has originally proposed a 
first-order BCMT because of the role of the long-ranged 
part of the Coulomb interaction. 1 ' However, the nu- 
merical result shows that the first-order transition takes 
place even with the onsite interaction only. In this cir- 
cumstance, it is desired to clarify basic properties of the 
BCMT and the FCMT in a unified way to further eluci- 
date the contrast. 

The Hubbard model has been studied for a long time 
as a minimal model to represent the essence of the Mott 
insulator and metals with their transitions. However, if 
the Hubbard model is defined on a bipartite lattice, the 
insulating gap is believed to open even at infinitesimally 
small onsite Coulomb interaction, thereby generic feature 
of BCMT cannot be studied in this case. 

When the lattice structure of the system is "geometri- 
cally frustrated" , meaning a nonbipartite structure with 
the appearance of frustration effects in the antiferromag- 
netic spin exchange, the BCMT is expected to occur at 
a nonzero onsite Coulomb interaction, which enables the 
study of BCMT and FCMT with their comparison on the 
same ground. The Hubbard model on the square lattice 
with nearest-neighbor (NN) and next-nearest-neighbor 
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(NNN) hoppings (HMSL) is a typical prototype for such 
a system, since this model contains the triangular net- 
work of hopping integrals, where the magnetic frustra- 
tion between the NN and NNN exchange interactions 
arises. 10 - ) This model was also proposed to be a minimal 
model of the high-T c superconductors in a certain param- 
eter regime near the Mott insulator. The LDA calcula- 
tions support that the cuprates are approximately de- 
scribed by the HMSL. 11, 12 ) Among the high-T c cuprate 
superconductors, (La,Sr) 2 Cu04 compounds are fitted by 
the parameters of the NNN hopping in unit of the NN 
hopping, —0.2, while the YBa2Cu307-<5 as well as Bi 
compounds are fitted by —0.3 <~ —0.5. On more general 
grounds, the nature of the electronic states around the 
Mott transition of the HMSL with nonzero NNN hop- 
pings has been a intriguing issue by its own right. In- 
deed, with the frustration effects, the ground state can 
be highly nontrivial, since spin entropy is not easily re- 
leased because of the competing magnetic interactions. 

Numerical methods offer potential tools for studying 
strongly correlated systems such as the Hubbard model. 
However, the phase diagram of the correlated metals 
and the Mott insulators in the Hubbard-type models has 
never been fully clarified until recently, partially because 
the geometrical frustration effect causes serious difficul- 
ties in numerical approaches. Recently, a new numer- 
ical algorithm, the path-integral renormalization group 
(PIRG), has been developed, which enables us to ob- 
tain the ground state accurately even in the system with 
strong frustration effects. 13 ' 14 ) A remarkable point is 
that the PIRG does not suffer from the negative-sign 
problem and can be applied to any type of Hamiltonian 
in any dimension. This method improves the ground 
state from the Hartree-Fock state by increasing non- 
orthogonal basis functions and performing the renormal- 
ization in the imaginary-time direction so that electron- 
correlation effects and quantum fluctuations are taken 
into account systematically. To reach the exact ground 
state of finite systems in a controlled way, the zero vari- 
ance limit is taken, which corresponds to the extrapola- 
tion of the truncated Hilbert space to the exact one. 

In this paper, we extend the algorithm of the PIRG 
to the grand-canonical ensemble. This grand-canonical 
path-integral renormalization group (GPIRG) method 
allows us to study the FCMT and BCMT within a single 
numerical approach in an efficient way. We determine the 
ground-state phase diagram of the HMSL in the plane of 
the chemical potential and bandwidth. A V-shaped Mott 
insulating phase is identified. We derive a general equa- 
tion, in which the slope of the metal-insulator transition 
line in the phase diagram is expressed by the ratio of 
the jump of the double occupancy to that of the electron 
density at the first-order transition point. We also de- 
rive an equation which connects the slope to the charge 
compressibility when the continuous transition occurs. 
By analyzing these relations and the chemical-potential 
dependence on the electron density, we reach a consis- 
tent picture for the V-shaped phase boundary and a 
sharp contrast of the FCMT with the first-order BCMT. 
The GPIRG also allows us to estimate the Coulomb- 
interaction dependence of the charge gap accurately. 



The organization of this paper is as follows: In §2, the 
framework of the GPIRG method is presented and tech- 
nical details for the implementation are described. In 
§3, the algorithm of the GPIRG method is examined by 
comparing with other methods such as the exact diago- 
nalization and the QMC methods. We show that the 
GPIRG is useful in calculating the chemical-potential 
dependence of physical quantities. The results of the 
HMSL by the GPIRG are reported in §4. The nature of 
the FCMT and the BCMT is discussed in detail. The 
summary is given in §5. 

§2. Grand-Canonical Path Integral Renormal- 
ization Group Algorithm 

In this section, we propose the grand-canonical path- 
integral renormalization group (GPIRG) method. The 
GPIRG is useful for studying the chemical-potential 
dependence of physical quantities as described below. 
Though this algorithm can be applied to a general form 
of electronic Hamiltonian on a lattice, we show the for- 
malism for the Hubbard model as a typical example. 
First we describe the key elements of the GPIRG method 
from §2.1 to §2.4: We introduce the particle-hole trans- 
formation in §2.1 and the path- integral operation in the 
Slater-determinant representation is explained in §2.2. 
The truncation and optimization procedure, and the 
extrapolation procedure of physical quantities are ex- 
plained in §2.3 and §2.4, respectively. The whole pro- 
cedure of the GPIRG is summarized in §2.5. Technical 
details for the implementation of the GPIRG are noted 
in §2.6 and §2.7. 

2.1 Model and Particle-hole Transformation 
The Hubbard model which we consider is 

where Ci„ {c icr ) is the annihilation (creation) operator on 
the i-th site with spin a and n ia = c\ a c ia in the A^-lattice 
system. The transfer integral is taken as tij = t for the 
nearest-neighbor sites and tij = t' for the next-nearest- 
neighbor sites. Throughout this paper, we take t as the 
energy unit. The total electron number is 

N 

2(7 

where (...) represents the ground-state expectation value. 
The canonical transformation 15 ) is introduced as 

Cfet — > Cfe, 

c-kl - 4- ( 2 ) 

In terms of the new operators, the Hamiltonian is rep- 
resented as 

H = H K + Hu- HL+^N, (3) 
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h k = ~Y^ ^3 (4 c j + c ] c i) + ( j - m) E 4 



{iJ) 



U 



The total number of electrons is related to the differ- 
ence of c and d particles as 



N 



N e = N + Y,(4 c i-4 d i)- 



(4) 



i=i 



The magnetization is given by 

N 

i=i 

= -N + M, 

where M is the total number of c and d particles, 

JV 



(5) 



M = Y,{4c i + d\d i ). 



i=l 

2.2 Projection procedure of kinetic and interaction 

terms in the Slater determinant basis 
2.2.1 path-integral operation 

In the GPIRG, numerical renormalization is performed 
in the imaginary-time direction. The ground state \ip g ) 
is obtained by 



\ip g ) = lim exp[-TH]\ip ), 



(6) 



where IV'o) is an initial state which is not orthogonal 
to \4> g ). Practically, following Feynmann's path-integral 
formalism, the projection procedure is performed by tak- 
ing sufficiently small A r : 

\ip g ) = lim lim (cxp[-A T iJ]) n |V>o), 

A T — iOb- >oo 

where we take r = A T n sufficiently large. 

The projection operator can be divided into the kinetic 
term and the interaction term approximately: 

e -A T (H K + Hu ) = e -A T H Ke -A T H v + q ^ (y) 

In this paper, we assume that the basis states are the 
Slater determinants. By using eq. (7) and the basis set 
of the Slater determinants, the projection to the ground 
state corresponding to eq. (6) is performed as described 
below. 

2.2.2 Slater determinant basis 

We use the notation \<p) to represent a Slater determi- 
nant. The Slater determinant which is a single-particle 
state is represented by the 2N x M matrix, <j>: 

M I IN \ 

w=n i°>> ( g ) 
k=i \i=i ) 

where c, = Cj for i = 1,...,N site and = g^-tv for 
i = N + 1, ...,27V site. The difference from the PIRG 
is that there appear the off-diagonal matrix elements 



[4>] ik for i = 1,...,N and k = N + 1,...,2N and for 
i = N + 1,...,2N and k = 1,...,N, which hybridize c 
and d particles in the GPIRG. 

We note that the off-diagonal elements of \(j)]ik cor- 
respond to the superconducting order parameter, i.e., 
{c\dk) — ( c i| c Lfe|) m the ground-canonical ensemble. 
As will be explained below, the GPIRG offers a system- 
atic improvement of a mean-field solution such as the 
Hartree-Fock or the variational Monte Carlo method. 
For example, we can set the BCS wavefunction as an 
initial state of the GPIRG and the fluctuation beyond 
it can be taken into account as the GPIRG calculation 
proceeds. In Appendix A, the BCS wavefunction with 
various symmetries is given in the GPIRG framework 
and we make some comments on it. 

We note that if we take the 2N x N Slater matrix, it 
is nothing but to specify the subspace with S z = from 
eq. (5). Hereafter we consider the S z = subspace and 
hence take the 2N x N Slater matrix, [<j>]ik- The con- 
straint S z = docs not restrict our calculation and even 
the ferromagnetic state can be represented in principle 
in the form of the magnetization in the xy plane. 
2.2.3 interaction-term projection 

The projection to the ground state is performed by 
using eq. (7) in the Slater-determinant basis. Though the 
interaction introduces a many-body term, the projection 
by the interaction term can be transformed into the sum 
of two Slater determinants by using the Stratonovich- 
Hubbard (SH) transformation. As seen in eq. (4), the 
total electron number is represented by the difference 
between c and d particle numbers. Then, we see that if 
the numbers of c and d particles change according to ji, 
the algorithm allows the grand-canonical framework at 
a fixed n. Hence, we introduce the SH transformation 
which hybridizes c and d particles as follows: 



exp 



-A T U (c\ Cl + d\d^j - c\ Cl d\d 
= ^ E cxp ^ s { c i d * + 4 C i) 



s=±l 



(9) 



where (3 = cos" 1 [cxp(-A T C//2)] for A T U > 0. Here, 
s = ±1 are the SH variables. A proof of eq. (9) is given 
in Appendix B. The above projection has the off-diagonal 
form with respect to c and d operators. 

We also note that the SH transformation with the di- 
agonal form for the attractive Hubbard model 16 ) is given 
as 

1 x > / „ A T U' 
exp 



exp 



A T Uc} Ci d]d 



E 

=±i 



-2as 



x exp 



where 



2as + 



A T U 



exp 



2as 



A T U 



dldi 



(10) 




for A T U > 0. 
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By using the SH transformation, the projection for the 
interaction term is performed as 

1 



expressed as 



exp 



ArUctadUi 



\<P) = , (\4> + ) + l<T» , (ii) 



where the sign in the right-hand side of eq. (11) specifies 
each SH variable s = 1 or — 1. In this way, the projection 
by the local interaction generates a sum of two Slater 
determinants. 

By the interaction-term projection on the i-th site, the 
right-hand side of eq. (11) is expressed as follows: 



M 



= n 



27V 



^ [I + X{+l)]..,[<t>] jlk c] 



k=l 



M 



2N 



10), 



10), 



where / is the 2N x 2N unit matrix and X(s) is the 
2N x 2N matrix. When the offdiagonal-type projection 
by eq. (9) is used in eq. (11), X(s) is given by 



[X(s)l 

[X(*)]N+i,i 



cos(/3s) — 1 
i sin(/3s) 



i X ( s )]N+t,N +I 



i sin(/3s), 

cos(/3s) — 1, 



and otherwise. Here, the local interaction projection is 
performed on the i-th site as in eq. (11). 

If we use the diagonal-type projection by eq. (10), X(s) 
is given by 



[X{s)] n , = cxp 



[X(s)} n , = cxp 



2as + 



2as 



A T U 
2 

A T U 



- l.for j = j' = i, 

- l.for j = j' = N + i, 



and otherwise. 

2.2.4 kinetic-term projection 

The kinetic-term projection is performed by multiply- 
ing the kinetic term of eq. (7) to the Slater determinant: 



cxpi-A^GMM = W). 



(12) 



Then, the projection of the single-body operator trans- 
forms a Slater determinant to another single Slater deter- 
minant. Here, fi p d is an artificially introduced chemical 
potential which is different from the real \i. We call /z p d, 
the pseudo chemical potential. By operating eq. (12) 
with several choices of ^ p d, the candidates for the ground 
state for given /z can be generated. Using /i p d instead of 
the real \i in the projection process helps escaping from 
a local minima with a wrong electron number separated 
by a potential barrier from the real ground state. 

2. 3 Truncation of the Hilbert space 

As seen in eq. (11) after the projection of the iV-sites 
interaction, exp[— A T Hu], an original single Slater deter- 
minant expands to the sum over 2 N Slater determinants. 
Then it is necessary to truncate the expanded Hilbert 
space to the optimal state for the ground state, which is 



= ^2w a \4> a ), 



(13) 



where L is the number of the optimized Slater determi- 
nants and w a is the optimized coefficient. To perform 
the truncation and optimization of the states, we select 
the set of the Slater determinants with the lowest en- 
ergy after each local-interaction projection, eq. (11) and 
the kinetic-term projection, eq. (12). This procedure is 
performed on the basis of the variational principle. 13, 14 - ) 
Namely, after each projection we solve the generalized 
eigenvalue problem 

L L 

Y,{<f> a \H\4> b }w b = E^{4> a \<t> b )w b , (14) 



6=1 



6=1 



and obtain the optimized set of ui b . To calculate 
(4> a \H\4> b ), the Wick's theorem can be used, since the 
basis set is constructed by the sets of the single-particle 
Slater determinants. 

2.4 Variance extrapolation 

If we can store the dimension L equal to the whole 
Hilbert space, eq. (13) gives the exact ground state. 
However, in practice the whole Hilbert space is not 
tractable for large system sizes. Hence, we use the ex- 
trapolation procedure by the energy variance. 13, 14, 17 -* 
This is based on the fact that the deviation from the 
ground state energy is proportional to the energy vari- 
ance 



Se oc Af 



(15) 



when the optimized state of eq. (13) is a good approxi- 
mation of the ground state. Here, Se is given by 

S E = (H) (H) g , 

where (H) and (H) g denote the ground-state energy with 
restricted and whole Hilbert spaces, respectively. The 
energy variance is defined by 



A; 



((g 2 )-(g) 2 ) 



By using the relation of eq. (15), the ground-state energy 
is obtained by the linear extrapolation to the A^ — > 
limit. The expectation value of the physical quantities 
such as equal-time correlation functions is also obtained 
by taking the linear extrapolation to the Ae — » limit. 

2.5 Whole procedure of the GPIRG method 

Until the previous section, we explained the key ele- 
ments of the algorithm. Here the whole procedure of the 
GPIRG is summarized. 

1) Let us start with the state, — X)a=i ^al^a)- 
When L — 1, the lowest-energy state = \(pi) is 

given by the Hartree-Fock solution. We can also set a 

variational basis state explicitly as an initial state of the 

GPIRG as shown in Appendix A. 

First we choose a basis state \<p a ) from L stored ba- 
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sis states {|0i), \ <h)> •••> I0i)}> which will be operated by 
exp[-A T H]. 

2) The projection by the kinetic term is performed as 

W=e- 4 ' fl « [ W|W, (16) 

where ^ p d is the pseudo chemical potential. In prac- 
tice, we set several ^ p d's around and the kinetic-term 
projections are performed by changing [i p d- After each 
operation of eq. (16), we solve the generalized eigenvalue 
problem in eq. (14) for the two basis sets: 

|0 2 ), |0 O _!), |0 a ), |0 O+1 ), |0 r ), (17) 

l^),!^),...,!^.!),!^),!^!),...,!^). (is) 

We always calculate the expectation value of the Hamil- 
tonian in eq. (14) by using fi, but not /i p d- Then, we se- 
lect either of the basis sets of eq. (17) or eq. (18), which 
gives the lower ground-state energy. By employing the 
pseudo chemical potentials, the change of the electron 
number is promoted toward the ground state for given /i 
after the GPIRG procedures of 2) and 3) described be- 
low. The kinetic-term projection by eq. (16) is done for 
each \<t> a ) for a = 1, . . . , L. 

3) The projection by the local interaction term is per- 
formed as 



(19) 



To make the electron number change toward the ground 
state for given fi, the off-diagonal-type projection in 
eq. (9) is employed. We solve the generalized eigenvalue 
problem in eq. (14) for the three basis sets: 

|01>, 102), |0a-l), \<f>a), |0a+l), -, 104 , (20) 
|0!), I^),...,!^.!),!^),!^!),..., \cj> L ), (21) 

|0 1 ),|0 2 ),...,|0 a _ 1 ),|0-),|0 a+1 ),...,|0 L ). (22) 

Then, we select the basis set among eq. (20), eq. (21), 
and eq. (22), which gives the lowest ground-state en- 
ergy. Taking either of |0+) or |0~) leads to the change 
of the total-electron number. The local-interaction-term 
projection is performed on the i = 1, . . . , N site. This 
procedure is done for each |0 a ) for a = 1, . . . , L. 

4) To increase the Hilbert space dimension, L, a new 
state is generated. Practically, we multiply the interac- 
tion term e ATUc i° id i di to a state \<f> a ) and select either 
state of |0+) or |0~). Here, off-diagonal-type projection 
eq. (9) is employed, since the hybridization between c and 
d particles promotes the change of the electron number 
in the grand-canonical ensemble. The state |0 a ) and the 
local site i to which eq. (9) is operated are selected ar- 
bitrary by using random number, respectively. We also 
select either state, |0+) or |0~) by using the random 
number. 

5) The procedure from 2) to 4) is repeated until the suf- 
ficient number of basis L is reached. Typically, it is suf- 
ficient to take L = 200 <~ 400. 

6) To reach the true ground state in finite systems, we ex- 



trapolate the physical quantiles such as the ground-state 
energy and correlation functions to the zero-variance 
limit as explained in §2.4. Namely, in the plane of the 
physical quantity (A) and the energy variance A^, the 
linear dependence appears in the case of the well con- 
verged states. Then, the extrapolated value by linear 
extrapolation is obtained finally. 

Here we note that in the GPIRG the error bar can 
appear by the two procedures: One is the variance ex- 
trapolation in each system size. If the linearity in the fit- 
ting as a function of the energy variance becomes worse, 
the error bars increase. The other comes from the ex- 
trapolation to the thermodynamic limit by assuming the 
finite-size corrections. The latter is a common origin of 
error bars to all the other numerical methods for finite- 
size systems, while the former is specific to this GPIRG 
method. The error bars in our analyses are given from 
the combination of these two types of error bars. 

In the interaction-term projection, we employ the 
diagonal-type projection in eq. (10) on the sites selected 
by the random number, in addition to the off-diagonal- 
type projection in eq. (9), both of which with an opti- 
mized combination make the convergence faster empir- 
ically. This is due to the fact that the diagonal-type 
projection gives the smaller variance by improving the 
basis state within the fixed electron number, while the 
off-diagonal-type projection tends to change the states 
with different electron numbers, which leads to the larger 
variance. 

We also note that in the GPIRG the negative-sign 
problem which often becomes serious in the QMC 
method does not appear, since the variational wave func- 
tion is explicitly given, \ip) = Y^=i w a\4>a) ■ Hence, the 
GPIRG can be applied to any lattice structures such 
as one-, two-, and three-dimensional systems with any 
boundary conditions. 

2.6 Implementation of GPIRG 

For implementation of the GPIRG procedure, the 
techniques developed in the QMC method 6 -* are avail- 
able. The matrix elements for Hamiltonian in cq. (14) 
are calculated from the single-particle Green's function 
by using Wick's theorem, because the Slater determinant 
is a single-particle state. The inner product of two Slater 
determinants is given by 



h) =det (*[0 a ][0b] 



(23) 



where [0 a ] is the 2N x N Slater matrix in eq. (23). The 
single-particle Green's function is defined by 

where i and j specify the site index from 1 to 2N and Cj = 
Cj for i = 1, N site and Cj = rfj_jv for i = N+l, 2N 
site. The 2N x 2N matrix of the Green's function is 
calculated as 



2N 2N 



H^LLW^^'Wj';, (25) 

i'=l j'=l 
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where 



In the case that (<j> a \ is updated in eq. (19), 



[g ab ] 



The physical quantities such as equal-time correlation 
function can be calculated by using the Wick's theorem 
and the Green's function, eq. (24). ^From the computa- 
tional point of view, the calculations for the determinant 
in eq. (23) and the inverse matrix in eq. (25) cost com- 
putational time proportional to N. However, this can be 
saved by taking the following procedures. 

2. 7 Update of inner product and Green's function 

Let us consider the case that \<j>b) is updated by the 
interaction-term projection in eq. (19): 



nk 



blnk 



2N 

E 



[I + X(s)} nl [M lk - (26) 



Hereafter we consider the case that the projection is per- 
formed on the i-th site as in eq. (19) and write j = N + i, 
until the end of this section §2.7. 

The updated inner product is given by 

(<f> a \4>' b ) = D{<P a \4> b ), (27) 
where D = det(J + GX) = det(J + XG) is given by 



D = -F ij F ji + (l + F ii )(l + F jj ), 



with 



Fin — GaX{ n -\- Gij Xj n , 

Fjn GjiX{ n -\- GjjXj n . 
The updated Green's function 

_ (</>a|4cm|$,) 



G' 



ab 



(28) 



(29) 



(30) 



{4>aW b ) 

is represented in the notation with [ ab ] omitted, as 
follows: 



G', 



G n 



where 



+ 5 ni 

+ S nj 
1 



?jm) , (31) 



G nrn — G nm ^ \{G n iXa -\- G n jX ji} 
x {1 + Fjj)G im - {-iy+iF lo G om } 

x {(1 + F u )G jm - {-ly^FjiGim}] 



(32) 



By using eqs. (27), (31) and (32), the updated inner 
product and Green's function G' can be obtained by us- 
ing G. Namely, if the Green's function G before updating 
is known, the updated inner product and Green's func- 
tion G' can be obtained by eq. (27) and eq. (31), respec- 
tively without calculating the determinant in eq. (23) 
and the inverse matrix in eq. (25). 



2JV 



n=l 



the updated inner product {(j>' a \(j>b) is obtained in 
cqs. (27), (28) and (29) with X replaced by *X. The 
updated Green's function 



G' 



Wa\<t>b) 

is obtained in the form with [ ab ] omitted as follows 



(33) 



G', 



Gn 



GniXa -\- G n jX{j j S t 



+ 



(34) 



where G nm is calculated in eqs. (28), (29) and (32) with 
X replaced by t X. 

By using eqs. (31) and (34), we can calculate the ma- 
trix clement {4> a \H\4>b) and (4> a \<Pb) from the Green's 
function of the previous step after each local-interaction- 
term projection, eq. (19). 

We note that these expressions for update of the in- 
ner product and the Green's function can be used in 
the general case where the interaction term affects the 
i and j sites. Namely, these are not restricted to the 
GPIRG method, but can be used in the PIRG and QMC 
methods: For example, after multiplication of the Slater 
determinants and the interaction term such as the ex- 
change interaction and pair-hopping term in the multi- 
orbital system, eqs. (27), (31) and (34) are available for 
update. 

§3. Examination of the GPIRG method 

The accuracy of the GPIRG method is examined by 
comparing with results by other methods such as the 
exact diagonalization, the QMC, and the PIRG methods. 

3.1 Comparison of energy 

First we compare the ground-state energy between the 
GPIRG and other methods. In Fig. 1 the extrapolation 
of the ground-state energy by the GPIRG with [i = 0.0 
for t = 1.0, t' = 0.0 and U = 4.0 on the 6 x 6 HMSL 
under the periodic boundary condition is shown by filled 
triangles. The total electron number N e versus(vs.) vari- 
ance for the corresponding data is shown in Fig. 2 by 
filled triangles, and we see that this state is converged 
into A^ e = 36 after the variance extrapolation. Since we 
consider the S z = subspace in the GPIRG in this pa- 
per, it turns out that the state has the electron number 
(Nf,Ni) = (18, 18). The ground-state energy after vari- 
ance extrapolation is estimated at E = —66.94 ± 0.05. 
For comparison, we show the PIRG results (open circle) 
by setting (N^N^ = (18,18) in Fig. 1. The ground- 
state energy by the PIRG after variance extrapolation 
is estimated at E = —66.92 ± 0.04. 14 ) A more accu- 
rate estimate indicated E = —66.88 with the error less 
than 0.01. 18 ) The number of states up to L — 320 are 
shown in both the GPIRG and the PIRG methods. The 
ground-state energies by the GPIRG and PIRG meth- 
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ods are consistent within the error bars by the variance 
extrapolation. The QMC data, E = -66.96 ± 0.07 14) is 
also plotted as the cross. We see that the GPIRG data 
and the PIRG data are consistent with the QMC data. 



-50 



QMC £=-66.96 + 0.07 




0.015 



Fig. 1. Extrapolation of the ground-state energy by the GPIRG 
(solid triangle) with fi = and to the zero-variance limit in 
the HMSL for t = 1.0, U = 4.0 on the N = 6x6 lattice under 
the periodic boundary condition. The PIRG data for (iVj, Ny) = 
(18, 18) is represented by open circles. The QMC result is shown 
by cross on the (H) axis. 
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Fig. 2. Extrapolation of the total electron number N e by the 
GPIRG for fi = —1.50 (open squares) and for fi = 0.0 (solid 
triangles) in the HMSL for t = 1.0, U = 4.0 on the N = 6x6 
lattice under the periodic boundary condition. 



We have also checked the GPIRG with other methods 
in larger system sizes. For example, when the ground- 
state energy for the zero-variance limit of the GPIRG 
with \i = 0.0 is written by E and that of the PIRG 
for (N-\,N\) = (N/2, N/2) is written by E', the relative 
error defined by (E-E')/E is estimated at 2.26 x 1(T 4 for 
the N = 8 x 8 system and 1.98 x 1(T 4 for the N = 10 x 10 
system for t = 1.0, t' = 0.0 and U = 4.0. Since it has 
been confirmed in ref. 14 ' that the PIRG data agree with 
the QMC data for TV = 8 x 8 and 10 x 10 systems within 
the relative error less than 0.3%, these results indicate 
that the GPIRG also provides reliable data comparable 
to these methods even in large-sized systems. 

The comparison in the energy between the GPIRG and 
the other methods are summarized in Table. I. We see 



that the relative error is less than 0.3%. 



3. 2 Comparison of other physical quantities 

Next we compare the GPIRG method with the exact 

diagonalization method for correlation functions. 

The equal-time spin correlations in the momentum 

space is calculated from 



S(q) 



1 

3N 



N 



5> 



(35) 



where Si is the spin operator of the i-th site and each 
element of the spin is calculated from 



Sf = \ [st 
S! = I K T - 



2 V l 1 



Cil + c'^di 



s-) = ; 



n il) 



c !t c u 



4l c iT 



The expectation value of eq. (35) is calculated by using 
the Wick's theorem and the Green's function, eq. (24). 
Namely, after the particle- hole transformation, eq. (2), 
the expectation value of operators is decomposed by the 
Wick's theorem. The difference from the PIRG case 
is that the off-diagonal term for j and J, operators be- 
comes finite in the GPIRG framework; for example, in 
the representation after the particle-hole transformation, 
{cid\)(dic\) can have a nonzero value. We performed 
the GPIRG with n = -1.25 for t = 1.0, t' = 0.0 and 
U — 4.0 on the 6x2 HMSL under the periodic bound- 
ary condition. All the extrapolations of the equal-time 
spin correlations are shown in Figs. 3 and 4. The results 
after the extrapolation to the zero- variance limit is shown 
in Table. II. By a proper choice of the chemical poten- 
tial, the GPIRG result indicates the particle number con- 
verges to N e — 10 which is nothing but (iV-f, iVjJ = (5, 5) 
electrons. The result of the exact diagonalization for 
(Nf,Ni) = (5,5) electrons is shown for comparison in 
Table. II. 
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Fig. 3. Extrapolation of the equal-time spin correlations 
S(q x ,q y ) to zero energy variance of the GPIRG method with 
H = -1.25 in the HMSL for t = 1.0, U = 4.0 on the N = 6x2 
lattice under the periodic boundary condition. 
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system 


GPIRG 


exact diagonalization and Monte Carlo results 


relative error 




E=(H) 


E> = (H) 


\(E-E')/E\ 


6 X 2, 5 T, 5 I 


-25.70 ± 0.02 


-25.6952 


0.00022 


6 x 6, 13 T, 13 I 


-58.45 ± 0.09 


-58.32 ± 0.072 


0.0022 


6 x 6, 18 T, 18 I 


-66.94 ± 0.05 


-66.96±0.07 


0.00031 



Table I. Comparison between GPIRG and exact diagonalization and Monte Carlo results in the ground-state energy for t = 1.0 and 
U = 4.0 for various-size systems under the periodic boundary conditions in x and y directions. 




Fig. 4. Extrapolation of the equal-time spin correlations 
S(q x ,q y ) to zero energy variance of the GPIRG method with 
fj, = -1.25 in the HMSL for t = 1.0, U = 4.0 on the N = 6x2 
lattice under the periodic boundary condition. 



are shown in Figs. 5 and 6. The results after the extrap- 
olation to the zero-variance limit by the GPIRG method 
are shown in Table. III. 
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Fig. 5. Extrapolation of the momentum distribution functions 
n(q x ,q y ) to zero energy variance of the GPIRG method with 
H = -1.25 in the HMSL for t = 1.0, U = 4.0 on the N = 6x2 
lattice under the periodic boundary condition. 



(q*,qy) 


GPIRG 


exact diagonalization 


relative error 


(1,0) 


0.0497 


0.0488 


0.018 


(2,0) 


0.0462 


0.0458 


0.009 


(3,0) 


0.0466 


0.0457 


0.019 


(0,1) 


0.2765 


0.277 


0.002 


(1,1) 


0.2807 


0.281 


0.001 


(2,1) 


0.2858 


0.286 


0.001 


(3,1) 


0.2791 


0.279 


0.000 



Table II. Comparison between GPIRG and exact diagonalization 
(ED) in the momentum distribution S(q x ,q y ) for t = 1.0 and 
U = 4.0 on the N = 6 X 2 lattice for (7V T , ATjJ = (5, 5) electrons 
under the periodic boundary condition. 



0.08 




The momentum-distribution function is defined by 



n(q) 



1 N 

— Y 

i -3 



c lt c iT 



a iq-(R 4 -Rj) 



(36) 



where is the vector representing the place of the i-th 
site. The expectation value of eq. (36) is calculated by 
using the Green's function, eq. (24). A relatively large 
error is seen for n(q) at (q x -,Q y ) = (3,0). This is re- 
lated to the fact that the amplitude itself is already small 
(= 0.0185). All the extrapolations of the momentum- 
distribution functions for the same system as in Fig. 3 



Fig. 6. Extrapolation of the momentum distribution functions 
n(q x ,q y ) to zero energy variance of the GPIRG method with 
H = -1.25 in the HMSL for t = 1.0, U = 4.0 on the N = 6x2 
lattice under the periodic boundary condition. 



In all the equal-time spin correlations and the momen- 
tum distributions, except n(3, 0), the relative errors are 
less than 3%. 
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GPIRG 


exact diagonalization 


relative error 


(0,0) 


0.9681 


0.9681 


0.0000 


(1,0) 


0.9606 


0.9601 


0.0005 


(2,0) 


0.9291 


0.9281 


0.0011 


(3,0) 


0.0157 


0.0185 


0.1783 


(0,1) 


0.0697 


0.06806 


0.0235 


(1,1) 


0.0448 


0.04513 


0.0074 


(2,1) 


0.0275 


0.02787 


0.0135 


(3,1) 


0.0223 


0.02281 


0.0229 



Table III. Comparison between GPIRG and exact diagonaliza- 
tion (ED) in the momentum distribution n(q x ,q y ) for t = 1.0 
and U = 4.0 on the N = 6 X 2 lattice for ( jVj , ) = (5,5) 
electrons under the periodic boundary condition. 



3.3 Change of total- electron number 

Next, to clarify how the grand canonical ensemble 
works in our GPIRG algorithm, we examine the relation 
between the chemical potential and the total electron 
number. In Fig. 7, the chemical potential fj, vs. filling 
n = N e /N by the QMC data 7 -* are shown as open squares 
for t = 1.0, t' = 0.0 and U = 4.0 on the TV = 6 x 6 lattice. 
Here, the chemical potential \x is calculated by 

E(N e ) - E{N' e ) 



N P -NL 



(37) 



where E{N e ) is the ground-state energy and n — (N e + 
N' e )/(2N). The open circles with dashed line represent \i 
vs. n for £7 = 0.0. The step of the dashed line appears at 
the closed shell, N e = 26 in the non-interacting case. The 
open square pointed by the filled arrow is the chemical 
potential calculated by using eq. (37) from the QMC 
data for the canonical ensemble 7 ) at N e = 26 and N' e = 
28 and the open square pointed by the open arrow is 
obtained from the QMC data at N e = 24 and N e = 26. 
In the GPIRG, the total electron number N e is obtained 
as the output after the zero-variance extrapolation for 
given [i. To check whether the correct N e is obtained 
for the given /j, in the GPIRG, we start from a state 
at L = 1 with N e ~ 24 as an arbitrary initial state. 
If the input parameter fi is located between filled and 
open arrows in Fig. 7, to be consistent with the QMC 
data, the converged state by the GPIRG should have 
N e = 26. In Fig. 2, the GPIRG results arc shown by the 
open squares for /i = —1.50 and we see that the GPIRG 
data are actually converged into the ground state with 
N e = 26. Next we set /i = 0.0, which is larger than the 
/i indicated by the filled arrow in Fig. 7 and performed 
the GPIRG starting from the same L = 1 state as the 
above case at /i = —1.50. The results are shown by the 
filled triangles in Fig. 2, where the converged state has 
N e = 36. This also reproduces the correct result, since 
at /i = 0.0, half filling should be realized. Namely, in the 
present N = 6 x 6 system for t' = 0, half filling, N e = 36 
should be realized for /i = 0.0. 

In Fig. 8, we show N e vs. the variance normalized by 
A# at L = 1 for various input parameters, /i. Here, the 
same L = 1 state is used for each simulation. We see 
that at fi = —0.8, N e is converged to N e = 36, while 




Fig. 7. Chemical potential fi versus filling n calculated by the 
QMC 7 ' (open squares) and the GPIRG (solid diamonds) in the 
HMSL for t = 1.0, U = 4.0 on the N = 6x6 lattice under 
the periodic boundary condition. The open circles with dashed 
line is for U = 0.0. The filled arrow represents fi calculated by 
QMC 7 ' for canonical ensemble at N e = 26 and N e = 28 and 
open arrow represents fi calculated from 7V e = 24 and N e = 26 
(see text). 



at fj, = —0.9, it converges to N e = 26. The closed- 
shell electron number N e = 26 is stable until [i = —1.80. 
While at fi = —1.9 it seems to converge to N e = 24. 

The chemical potential vs. filling obtained from Fig. 2 
and Fig. 8 are plotted by solid diamonds in Fig. 7. We 
see that the GPIRG results are consistent with the QMC 
results. As seen from Fig. 7, the electron number sensi- 
tively and correctly converges with the resolution of 0.1 
for the input parameter fi in the GPIRG method. As the 
system size increases, the energy gaps of the shell struc- 
ture in finite-size systems of course become smaller. This 
makes the potential barrier small and the electron num- 
ber tends to change more continuously in the GPIRG 
calculation. On the other hand, in the canonical frame- 
work such as the PIRG, 13,14 ) the chemical potential is 
calculated from the subtraction of the ground-state en- 
ergies with different electron numbers in eq. (37). This 
gives larger error in the larger system sizes. Hence, the 
GPIRG is useful to calculate the chemical-potential de- 
pendence of the physical quantities such as the charge 
gap and the jjL-U phase diagram which will be shown in 
84. 



§4. Nature of Filling-Control and Bandwidth- 
Control Mott Transitions 

By using the GPIRG method we study the nature 
of the FCMT and the BCMT on the HMSL. We fo- 
cus on t'/t = -0.2 where the first-order BCMT oc- 
curs, 9 ) though the GPIRG can be applied to the larger- 
frustration regime. The GPIRG calculation is performed 
up to N = 10x10 lattices under the fully periodic bound- 
ary condition. 

4-1 Ground-state phase diagram in plane of chemical 
potential and band width 
By performing the GPIRG with the control parame- 
ters n and U, we obtain the total electron number A^. 
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Fig. 8. Extrapolation of the total electron number We by the 
GPIRG for fi = -0.8 (solid triangle) -0.90 (cross), -1.0 (open 
square), —1.70 (open triangle), —1.8 (open circle), —1.9 (solid 
circle), and —2.0 (open diamond) in the HMSL for t = 1.0, 
U = 4.0 on the N = 6x6 lattice under the periodic boundary 
condition. 



Then we can construct the ground-state phase diagram 
by plotting the boundary between the states with differ- 
ent electron numbers in the plane of fj, and U. Figure 9 
shows the ground-state phase diagram for t = 1.0 and 
t' = -0.2 on the N = 4 x 4 lattice. For U = the closed 
shells are realized at N e = 10, 14 and 22. The degener- 
acy of fi is lifted by switching on U. The solid line is the 
least-square fit of each boundary. Each area between the 
boundary lines corresponds to the state with each N e . 
The central triangle area is N e = 16, i.e., half filling. 
For U ^ 5 the boundaries between N e = 14 and 16, and 
N e = 16 and 18, are remarkably linear with U. This is 
expected in the large U regime, but it holds even close 
to the critical value U ~ 2.8. The linear opening of the 
gap basically comes from the energy cost to add an elec- 
tron (a hole) to half filling. The linear-fitting lines of the 
boundary between iV e = 14 and 16 and the boundary be- 
tween N e = 16 and 18 for U = 2.8, 3.0, 3.3, 3.5, 3.7, 4.0, 
4.5, 5.0, 6.0, 7.0 and 8.0 crosses at {p., U) = (-0.19, 2.80). 
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Fig. 9. Ground-state phase diagram in the plane of fi and U for 
t = 1.0 and t' = -0.2 on the the N = 4 X 4 lattice. The solid 
lines are the least-square fit of each boundary between states 
with different electron number. 




By performing the same procedure in larger-sized sys- 
tems, we construct the ground-state phase diagram in 
the plane of \i and U in Fig. 10. Here, each boundary 
between different N e states around n = 1 is shown for 
the system size of TV = 4 x 4 (black diamond) as well as 
for N — 6 x 6 (blue square). For TV = 8 x 8 (green circle) 
and N — 10 x 10 (pink triangle), the boundary between 
n = 1 and n / 1 closest to n = 1 is shown. The least 
square fit of each data is also plotted for the systems 
with N = 4 x 4 (solid line), 6x6 (blue line), 8x8 (green 
dashed line), and 10 x 10 (pink line). To estimate the 
boundary of the insulator phase for N — > oo, we extrap- 
olate the chemical potentials between the n = 1 state 
and the n ^ 1 state closest to n = 1 by using the fitting 
function 

M (iV) = (38) 



Fig. 10. Ground-state phase diagram in the plane of fi and U for 
t = 1.0 and t' = —0.2 on the square lattice. The data represents 
the boundary between different electron numbers in the N = 
4x4 (black diamond), 6x6 (blue square), 8x8 (green circle) 
and 10 X 10 (pink triangle) systems. Each line represents the 
least-square fit for each boundary in the JV = 4x4 (black line), 
6x6 (blue line), 8x8 (green dashed line) and 10 X 10 (pink line) 
systems. 



The particle number closest to half filling realized in 
Fig.10 is N e = 32, 62 and 98 for 6 x 6, 8 x 8 and 
10 x 10 lattices, respectively for the hole doping. For 
the electron doping, N e = 38, 66 and 102 in the same 
order as above. We employ the above finite-size scaling 
form because it well fits the data and partly because the 
Hartree-Fock gap equation also follows this form. 7 ) The 
least square fit by the form eq. (38) of \x for N — > oo at 
U = 4.0, 5.0, 6.0, 7.0 and 8.0 are shown by open circles in 
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Fig. 11. Ground-state phase diagram in the plane of fj, and U 
for t = 1.0 and t' = —0.2 on the square lattice in the thermo- 
dynamic limit. The solid lines represent the least-square fit of 
metal-insulator transition for U = 4.0, 5.0, 6.0, 7.0 and 8.0. Gray 
area represents the Mott-insulator phase in the thermodynamic 
limit. The open diamonds represent chemical potentials for half 
filling for U = 0.0 and U = 2.0 in the thermodynamic limit. By 
connecting those points to the BCMT point, the gray-dashed 
line represents the half-filling density in the metallic phase. The 
errors produced by the system-size extrapolation are not shown 
(see text). 



Fig. 11. 

The black lines in Fig. 11 represent a linear fitting 
of these open circles, which well represent the metal- 
insulator boundaries in the thermodynamic limit and 
two boundaries meet at (n,U) = (—0.38,3.15). The 
gray area inside the thick solid lines is the Mott insulator 
phase. Although the error bars arising from the system- 
size extrapolation become relatively large for small U 
(3.3 ^ U < 4) regime (see also Fig. 10), they are typically 
less than 0.1. In Fig. 11 we omit plotting the error bars 
arising from the system-size extrapolation. The open di- 
amonds represent chemical potentials for half filling for 
U = 0.0 and U = 2.0 in the thermodynamic limit. The 
gray-dashed line connects those points to the bottom of 
the Mott insulator phase, namely, the BCMT point. We 
see that the Mott insulator phase closes at U c = 3.3 in 
Fig. 11. This seems to be consistent with the PIRG re- 
sult at n = 1, in which a quite independent measurement, 
namely, the jump of the double occupancy indicating the 
first-order BCMT appears at U c = 3.25 ± 0.05. 9 > The 
metal-insulator boundary except the BCMT point is the 
boundary of the filling-control Mott-transition (FCMT). 
It should be noted that the phase boundary remarkably 
has a V-shaped structure rather than a U-shaped one. As 
in the black dashed curve in Fig. 11, the phase boundary 
of FCMT very close to the BCMT appears to show a 
small deviation from the linear fitting, this uncertainty 
coming from the uncertainty of the size extrapolation is 
small. It is not clear whether the phase boundary has 
a round shape (U-shape structure) near BCMT, or V- 



shaped until BCMT. However, the overall shape of the 
phase boundary is well represented by the V shape and 
the round structure is limited to a very narrow region if 
it exists. This remarkable V-shape structure will further 
be discussed in §4.3. 

4-2 Carrier dependence of chemical potential and 
charge compressibility 
Though we have obtained the slope SU /<5/i as in 
Fig. 11, the slope itself does not tell us the order of 
the metal-insulator transition. In the literature, the 
FCMT is identified as the continuous transition with a 
diverging compressibility by QMC studies, 7 ' 8 ' while the 
BCMT shows the first-order transition with a jump of 
the double occupancy, which has been obtained from the 
PIRG study. 9 ) Before discussing this contrast, we ap- 
ply GPIRG method to confirm whether the continuous 
character of the FCMT is stable and reliable, since the 
GPIRG method is an independent technique from the 
QMC method. To identify the order of the transition 
we calculate the filling dependence of the chemical po- 
tential. In Fig. 12(a) we show the chemical potential /i 
vs. filling n = N e /N for U = 4.0. The dashed line, 
the dash-dotted line, and the solid line represent the \x-n 
lines calculated by the GPIRG for the N = 4 x 4, 6 x 6 
and 8x8 systems, respectively. We plot the open sym- 
bols at the middle point of each step for the N — 4 x 4 
(open diamond), 6x6 (open square), and 8x8 (open cir- 
cle) systems. We also plot the data by the filled symbols 
calculated from eq. (37) by the PIRG in the closed-shell 
structure for the N = 4 x 4 (filled diamond), 6x6 (filled 
square), and 8x8 (filled circle) systems. As seen in 
Fig. 12(b) all symbols are fitted well by the function 



(39) 



with S = 1—n. The gray lines in Fig. 12(a) and Fig. 12(b) 
are the least square fit with fitting function, eq. (39), 
which gives a = -12.26 ± 0.27 and \i c = -0.76 ± 0.01. 
This indicates that the charge compressibility has the 
form, 



On 



n)~ 



(40) 



away from half filling as shown by Furukawa and 
Imada. 7 -* If the filling n becomes discontinuous as a func- 
tion of fiy the- first-order transition occurs. Namely, un- 
stable n leads to the inhomogeneous ground state which 
has the hole-rich and -poor regimes. On the contrary, 
the fitting of eq. (39) in Fig. 12 seems to show Xc — > °o 
for n — > 1 and Xc = at n = 1. In Fig. 12, however, the 
closest data point to half filling is n = 0.9375. Hence, 
we safely conclude that the phase separation does not 
occur for the carrier density larger than 6% for U = 4.0. 
Absence of the phase separation in the Hubbard model 
has been reported in refs. 7, 20 - ) for the no- frustration case, 
t' = as well as for a frustrated case at t' = O.2. 8 ) Our 
result is perfectly consistent with the QMC result and a 
and ii c show quantitative agreement with the result in 
Ref. 8 ) 

When the first-order transition occurs at the BCMT, 
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Fig. 12. (a) Chemical potential fi vs filling n for t = 1.0 and 
t' = —0.2 on the square lattice, (b) fi vs <5 2 . Dashed line, 
dash-dotted line, and solid line represent TV = 4 X 4, 6 X 6, and 
8x8 systems, respectively, obtained by GPIRG method. Open 
symbols are plotted at middle point of each step of fi-n line. Solid 
symbols arc obtained by PIRG method for closed shell structure. 
Gray curves in both panels are least square fits. 



one could expect that the phase separation also occurs at 
the FCMT in general. The numerical results discussed 
above, however, show rather contrasted behavior. The 
BCMT shows clear first-order transition, while FCMT is 
most likely continuous with a diverging compressibility 
at the transition point. 

In the framework of the commensurate AF mean- 
field approximation (Hartree-Fock approximation) (see 
cqs. (49) and (50)), the AF metal appears by the car- 
rier doping to the AF insulator and the phase separation 
between the paramagnetic metal and the AF metal oc- 
curs: For t' = —0.2 and U = 4.0, the phase-separated 
region appears for 0.69 < n < 0.85 in the hole-doped 
case. Namely, the phase separation appears for the in- 
terval of filling, 16%. The first-order BCMT occurs at 
t/F F = 2.064 at n = l 21 ) in the mean-field approxi- 
mation, 21 ~ 23 ) where the jump of the double occupancy 
is estimated at 0.04. By the PIRG calculation at n = 1 
the jump of the double occupancy is estimated at 0.04 at 
U c = 3.25±0.05. 9 ) Though the jump of the double occu- 
pancy at n = 1 by the PIRG method is comparable to the 
mean-field approximation, the region of the phase sepa- 
ration at the FCMT should be largely reduced from the 
mean-field result if it existed. This type of large phase- 
separation region in the mean-field approximation disap- 
pears in our result after considering the quantum fluctu- 
ation effect by the GPIRG. In addition, the phase sep- 
aration in the mean-field approximation occurs between 



a paramagnetic metal and antiferromagnetic metal. The 
FCMT itself is a continuous transition even in the mean- 
field solution, which is the same as our result, although 
the antiferromagnetic metal seen in the mean-field ap- 
proximation is likely to become absent in metals again 
because of the quantum fluctuations. 

We also comment about the contrast between the 
BCMT and the FCMT from viewpoint of the bound- 
state formation. Since the electrostatic attractive force 
works between a doublon (a site occupied by two elec- 
trons) and a holon (a site unoccupied by electrons), the 
bound state between a doublon and a holon is formed 
at the BCMT. The Mott insulating state is indeed re- 
garded as the phase where the doublons and holons all 
form bound states in pairs. Here we recall that in the 
dislocation- vector system 24, 25 ' and the Coulomb-gas sys- 
tem, 26 ) the Kosterlitz-Thouless transition changes to the 
first-order transition as the ratio of the binding energy 
of vortex cores and the core energy exceeds a thresh- 
old. The same mechanism may work in the BCMT and 
FCMT. Namely, in case of the BCMT, relatively large 
attractive force between a holon and a doublon induces 
the first-order phase transition. In the case of the FCMT, 
however, the mechanism is quite different. The transi- 
tion in this case is controlled by an interaction between 
two holons (or between two doublons). Apparently, the 
attraction between two holons should be much smaller 
than the attraction between a holon and a doublon. In 
the present system, if the attractive force works among 
holons (doublons) at the hole (electron)-doped Mott in- 
sulator, the bound state is formed between holons (dou- 
blons) and the first-order FCMT may occur. Our calcu- 
lation for U = 4.0 shows that the binding energy between 
holons is weak enough so that the transition becomes 
continuous or at most the discontinuity is very small. A 
sharp contrast between the FCMT and BCMT is natu- 
rally understood from this completely different origin of 
the attractions. 

4-3 Thermodynamic analysis 

To understand the contrast between the BCMT and 
the FCMT in detail, we derive a relation between the 
slope of the metal-insulator-transition boundary in the fi- 
ll phase diagram and some other thermodynamic quan- 
tities. Let us consider the expansion of the ground-state 
energy by fj, and U. 



E(fi + 5fi,U + SU) 
+0((Srf,(SU) 2 ). 



v ^ + \dU 



SU 



(41) 



In the /i-U phase diagram, along the metal-insulator- 
transition boundary, the energies of the metal and insu- 
lator phases may not change: 

E T {^U) - E M {H,U) 

= E I (j i + Sn, U + SU)- E M (fi + 5n,U + SU) 

= 0, (42) 
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where subscript I(M) represents that the expectation 
value is taken in the insulating (metallic) ground state. 
If the first-order metal-insulator transition occurs, the 
slope of the transition line in the \i-TJ phase diagram is 
determined by the ratio of the jump of filling and dou- 
ble occupancy. Namely, from cq. (41) and eq. (42) the 
following equation is derived: 

6U m - n M ^ 



where 
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At the BCMT point SU/Sn = if the first-order tran- 
sition occurs, since the numerator of eq. (43) becomes 
zero because nj—1 and um = 1 arc both satisfied, but 
the denominator is finite. 

In the case of the continuous metal-insulator transi- 
tion, let us consider the expansion of filling n: 



n(p + 5n,U + 5U) 
+0({8 i if 1 (W) 2 ) 
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By a parallel discussion to the above derivation, the slope 
of the metal-insulator boundary in the u-U phase dia- 
gram is derived as 

SU X (47) 



( dnM_\ ' 
V dU I p, 



where \ c is the charge compressibility 

' dn M " 

Xc = 



, , (48) 

in the metallic phase and {dni/dfi)u = = (dni/dU)^ 
is used. 

From the relations cq. (43) and eq. (47), it turns 
out that the V-shaped structure implies a quite differ- 
ent character between the bandwidth- and filling-control 
transitions. In fact, if one assumes the first-order tran- 
sition for the bandwidth control at the corner of the V 
shape as we indeed observed, dU/d/i is not well defined 
at the corner (see U c in Fig. 13(a)). At this corner, we 
have observed a jump in D, namely a nonzero Dj — Dm 
while the filling should not show a jump at this corner. 
This implies dU/du = 0, which does not contradict the 
character of the corner. 

However, if the first-order character is also retained 
in the edge of the phase boundary as drawn by thick 
lines in Fig. 13(b), namely for the filling-control transi- 
tion, the V shape together with eq. (43) results in an 
emergence of the abrupt jump in ni — um in the region 
above U c but infinitesimally close to the corner, which is 
unlikely. This is also understood in Fig. 13(c) which is 



the D-n phase diagram corresponding to Fig. 13(b). The 
abrupt jump in nj — tim at U infinitesimally close to U c 
in Fig. 13(b) leads to a nonzero, finite phase-separated 
region at U ~ U c in the D-n phase diagram as in the 
coexistence of R and Q (or R and P) in Fig. 13(c). At 
U = U c the ground-state energy at P, R and Q are de- 
generate also with that at S, because of the first-order 
nature of the BCMT. Namely, the double-well structure 
in the energy as a function of n (local minima at nj and 
Um) forces us to allow three different metal phases sep- 
arated by first-order transitions as shown as P, Q and 
S in Fig. 13(c). Then at U c , a metal with the density 
at P or Q is realized at the same chemical potential as 
the point R and S. Now the same density of metal as 
P and Q has to be realized also at A and B when the 
chemical potential is shifted from S, because the points 
A and B are realized from S by changing the chemical 
potential adiabatically. Then at the same U, the same 
density of metal (namely P and A ( as well as Q and B) 
) is realized at different two chemical potentials, which is 
unphysical. Therefore the case Fig. 13(b) is unlikely to 
be realized. The above discussion concludes that, the V- 
shaped metal-insulator transition lines together with the 
first-order BCMT is not compatible with the presence of 
the first-order FCMT near the BCMT. 

If the transition converts to the continuous character 
for U infinitesimally close to U c in Fig. 13(a), it implies 
a diverging compressibility and diverging dn/dll with a 
finite ratio for FCMT. This is certainly possible. There- 
fore, the V-shaped structure is consistent with the abrupt 
conversion of the first-order character observed in the 
bandwidth-control transition to the continuous charac- 
ter of the filling-control transition supported from the 
data in the previous subsection. 

In the case of the U-shaped structure, it is possible 
that the first-order transition occurs at the bottom U c 
as in Fig. 13(d), and also that the first-order transi- 
tion is retained until the critical points U' c (> U c ) as 
drawn by thick lines in Fig. 13(e). The latter is because 
dU / du changes from zero to finite values continuously as 
U changes from U c . Namely, as U changes from U c to U' c , 
Di—Dm changes from a finite value to zero, and Ui—Hm 
changes from zero to finite values and finally becomes 
zero, continuously. In the corresponding D-n phase dia- 
gram, Fig. 13(f) with the same notations as Fig. 13(c), 
we see that there exists no unphysical phase separation 
within the metallic phase in contrast with Fig. 13(c). In 
this case, the phase-separated region terminates continu- 
ously at both ends of U' c and U c , and the insulator phase 
is realized for U > U c . This argument does not in prin- 
ciple exclude the possibility that U' c becomes infinity. 

In our calculation, within the numerical accuracy, we 
do not exclude the possibility that the V shape could 
in reality have a U shape corner structure in the very 
vicinity of the BCMT with dU /du being well defined ev- 
erywhere. In this case, we do not need to assume an un- 
physical jump for the filling control as mentioned above. 
Even though, a rather sharp change observed at the "cor- 
ner" of the phase boundary supports a rapid change of 
the character between two routes of the transitions. 
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Fig. 13. Upper panel: Possible schematic phase diagrams with a 
V-shaped metal-insulator boundary, where the first order transi- 
tions are indicated by solid circles or thick lines while thin lines 
represent continuous transitions, (a) shows the first-order tran- 
sition only at the corner while continuous in the edges , and (b) 
shows the case with the first order transition both at the corner 
and edges (in the part drawn by thick lines), (c) —D-n phase 
diagram corresponding to the case (b). Lower panel: Possible 
schematic phase diagrams of U-shaped metal-insulator boundary 
with the same notations as the V-shaped cases, (d) shows the 
case with the first-order transition only at the bottom, and (e) 
shows the case with the first-order transition both at the bottom 
and edges (in the part drawn by thick lines), (f) —D-n phase 
diagram corresponding to (e). In (c) and (f), the insulator phase 
is represented by thick line at n = 1, and the metallic phase is 
represented by the gray area. — D = 0.25 is taken along the n 
axis. The miscibility gap due to the phase separation is given 
by the white area with dashed curves, which connect the coex- 
isting metallic and insulating points. Namely, the two ends of a 
dashed curve represent the jump of D and n at the first-order 
transitions for a given U < U' c . U c denotes the BCMT point 
and U' c denotes the possible critical point where the first-order 
transition converts to the continuous transition. The case (b) 
(and (c)) is unlikely to be realized. The present numerical result 
supports the case (a) (see the text). 



We are now left with the possibilities of (a), (d), and 
(e) in Fig. 13, where the round part is limited to the 
region close to the corner if (d) and (e) apply, because 
the overall V shape is clear. In fact, the analysis of the 
charge gap in the next subsection more precisely shows 
that the gap opens linearly with U — U c . Therefore the 
phase diagram is most likely to belong to the class (a) 
in Fig. 13, the FCMT is continuous while BCMT is the 
first-order transition. We also analyze whether the first- 
order transition can be consistent with the slope in the 
V-shaped structure away from the corner. The slope 
SU/S/u, is roughly 2.4 in the hole-doped part. When the 
first-order transition occurs, and the jump of the double 
occupancy has a similar value to the jump at BCMT (~ 
0.04), the phase separation should extend to the doping 
concentration ~ 0.1. However, this is clearly not the case 
from our numerical results, and the phase separation is 
limited at most to the region less than 0.06. Rather 
sharp difference in the character of FCMT from BCMT 
becomes recognizable also from such analyses. 

In the T-shaped structure case, such that the metal- 
insulator transition line behaves as \i ~ (U — Uc) 7 with 



7 > 1, the situation is the same as the V-shaped case. 
The essential singular form such as /i ~ exp(— 2itJ t/U) 
is known to be realized for t' = in the HMSL 19 ) with 
the metal- insulator transition at U = 0. We arc unaware 
of the type of transitions in a concrete model except the 
case U c = 0. If the system with this form with the first- 
order BCMT point U c > exists, it is classified to the 
same class as discussed above. Therefore, the first-order 
transition is likely to occur only at the corner. 



4-4 Charge excitation gap 
The charge gap is defined by 
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where fj,((2N' - 1)/N) = {E(N',N') - E(N' - 1,7V' - 
l)}/2 and E(N^,N^) is the ground-state energy with 
number of up(down) spin, iVj(iVt). 

Figure 14 shows the charge gap calculated by the 
GPIRG method. Solid squares represent the charge 
gap, A c in the thermodynamic limit. Note that, in the 
GPIRG method, the charge gap is determined by the dis- 
tance in jU between two phase boundaries in Fig. 11. For 
the extrapolation, we have used the scaling function 

A c (A0 = A c + -^=, 

as in the Hartrce-Fock AF gap equation. 7 ) Since the 
opening of the gap A c is very well described by a linear 
U dependence, the solid line is fitted by the linear line of 
A C (N -> oo) for U = 3.5, 4.0, 5.0, 6.0, 7.0 and 8.0. From 
this fit, the critical value U c is estimated at U c = 3.23. 
In Fig. 14, the gray diamond represents U c — 3.25 ±0.25 
at which the jump of the double occupancy appears at 
n = 1 in the PIRG calculation. 9 ) These two independent 
estimates agree each other excellently. The charge gap 
grows from continuously even at the first-order metal- 
insulator transition point U c and shows remarkably a 
linear dependence on U. This continuous growth of the 
charge gap from zero is rather obvious because the first- 
order transition indicates a level crossing of the metallic 
and insulating states, and the energy difference between 
the metallic and insulating states increase continuously 
from zero after the level crossing, while the charge gap in 
the insulating side must be smaller than this difference. 

In the Hartree-Fock approximation, 21 ~ 23 ) the AF 
bands are given by 

£ - ± (k) = i[ £ (k) + £ (k') 



± { £ (k)- £ (k')}Wi + 



U 2 m 2 



{e(k)- £ (k')} 2 



where (n ia ) = (n + ct(-1)I*Ito)/2 and k' = k + Q with 
Q = (it, 7r). The self-consistency equations are 
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where the summation for k is taken in the lst-Brillouin 
zone and /(e) is the Fermi distribution function. For 
t' = —0.2 the indirect gap appears between the AF 
bands. Namely, the charge gap A^ F is calculated from 
the energy difference between the bottom of the upper 
band e~(7r/2, n/2) and the top of the lower band e + (7r, 0) 
by using the solution of eqs. (49) and (50), \i and to, for 
n = 1: 

The charge gap A^ F is shown by the dashed line in 
Fig. 14. Here, the order parameter shows the jump 
from to = to 0.391 at the first-order-transition point 
/7,P F =2.064, 21 ' but the charge gap appears from to fi- 
nite values continuously for U > U^ F . Comparing A^ F 
with A c obtained by the GPIRG method, we see that U c 
is larger than U^ F and A c is reduced from A^ F by the 
electron-correlation effect: for example, A c /Aj? F = 0.32 
at U = 4.0 and 0.58 at U = 8.0. 

It is interesting to estimate the parameters of the 
HMSL as the effective model for the cuprates: Optical 
conductivity data of La2Cu04 27 ' indicate the charge gap 
2A C of the order of 1.5 - 2 cV. If the HMSL describes 
low energy physics of the material, the NN hopping is 
estimated at t — 0.43 eV according to ref. 11,12 ' Hence, 
the value of U/t is estimated at 7.8 ~ 9.3 from Fig. 14. 
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Fig. 14. (/-dependence of the charge gap for t = 1.0 and t' = 
—0.2 in N = 4 X 4 (open circle), 6x6 (open diamond), 8x8 
(open square), and 10 X 10 (open triangle) systems. Solid squares 
represent the charge gap for N — ► oo. Solid line is least square 
fit of A C (N — » oo) (see text). Gray diamond is located at U c = 
3.25 ± 0.05 determined by the jump of the double occupancy at 
n = 1 by PIRG. 9 - 1 Dashed line represents the charge gap A^ F 
by the AF mean-field approximation. 



Since the exponent of the charge gap a defined by 
A c oc \{U — U c )/t] a is consistent with a = 1 both in the 
Hubbard model in two dimensions and in the Hartree- 
Fock result, it may be a universal feature irrespective of 
the dimensionality, except a special case of U c = in the 
perfectly nested lattice. This universality of the expo- 
nent is intriguing because the linear opening of the gap 
also means the V-shaped phase boundary and supports 
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the first-order transition for the BCMT coexisting with 
the continuous transition for the FCMT. The opening of 
the Mott insulating gap as a function of (U — U c )/t may 
be compared with the experimental results in the litera- 
ture. In the perovskite compounds, Ri-^CaajTiOa, the 
opening of the charge gap has been studied by Katsufuji, 
Okimoto and Tokura. 2 ' 28 ' Here R represents rare earth 
ions to control effective (U — U c )/t. The optical gap 
as well as the activation energy in the resistivity show 
an opening of the gap well scaled by a linear function 
of (U — U c )/t. Our calculated result is quite consistent 
with this observation. Although a narrow doping region 
of the doped system in real three-dimensional perovskite 
may be under influence of impurity localization and also 
of the antiferromagnetic-metal phase, it turns out that 
the linear opening of the charge gap is tightly connected 
with a singularly different characters between FCMT and 
BCMT. When the critical point of the Mott transition 
characterized by the compressibility divergence occurs 
close to zero temperature as supported from the contin- 
uous character of FCMT, the character of the Mott tran- 
sition has to be analyzed not as the classical Ising-type 
transition 29 ^ but as the quantum phase transition. If the 
critical point is normally characterized by the vanishing 
k 2 dispersion replaced with the dominant fc 4 dispersion, 
this dynamical exponent z = 4 leads to the compress- 
ibility as k oc 8 1 ^ z / d with the doping concentration 5 in 
d- spatial dimensions in the scaling analysis. 30 ' In two 
dimensions, this leads to k oc S^ 1 , which is consistent 
with the numerical as well as experimental results. 2 ' In 
three dimensions, this leads to n oc <5~ 1//3 . The diver- 
gence is much weaker than the two-dimensional case and 
rather difficult to see the compressibility divergence ex- 
perimentally as compared to two-dimensional systems. 
This is consistent again with the experimental indica- 
tions, 31 ' although the surface effect in the photoemission 
experiments has to be carefully considered. 

§5. Summary 

We have studied the ground-state properties of the 
Hubbard model on the square lattice (HMSL) with 
nearest-neighbor and next-nearest-neighbor hoppings. 
The computation by our new algorithm has made it pos- 
sible to draw a phase diagram of the Hubbard model in 
the parameter space of the interaction U, chemical po- 
tential (i and the band structure (or amplitude of frus- 
tration). 

Our new method is summarized as follows: To study 
the bandwidth-control Mott transition and the filling- 
control-Mott transition by a unified approach, we have 
developed an algorithm of the grand-canonical path- 
integral renormalization group (GPIRG) method. To 
treat the system in the grand canonical ensemble, the 
particle-hole transformation is applied to the Hamilto- 
nian and the basis states. To reach the ground state 
for the chemical potential /z, the Stratonovich-Hubbard 
transformation which hybridizes up and down electrons 
is introduced. To avoid the trapping in a local minimum 
with a specific electron number, it is efficient to use the 
pseudo-chemical potentials (fictitious chemical potential) 
in the kinetic-term projection, which arc different from 
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the original chemical potential /U. This algorithm does 
not suffer from the negative-sign problem, which often 
becomes serious in the quantum Monte Carlo method 
and can be applied to any Hamiltonian in any lattice 
structure and dimension with arbitrary boundary condi- 
tions. The GPIRG results in the ground-state energy and 
the physical quantities such as the equal-time spin cor- 
relation and the momentum distribution function show 
good agreement with the results by the exact diagonal- 
ization and the quantum Monte Carlo methods. The 
GPIRG is shown to be useful to calculate the chemical- 
potential dependence of physical quantities. 

Our calculated results by the GPIRG for the two- 
dimensional Hubbard model is summarized as follows: 
The ground-state phase diagram of the HMSL in the 
plane of the chemical potential n and the interaction U 
is precisely determined. The contrast between the band- 
width and filling-control Mott transitions clarified in the 
literature are more firmly confirmed on the basis of a uni- 
fied approach here, in the following fashion: The carrier- 
density dependence on the chemical potential indicates 
that the phase separation does not occur (for example, 
at least for carrier doping larger than 6% at U = 4.0), 
implying the continuous character of the transition with 
enhanced fluctuations near the transition as seen in the 
diverging charge compressibility. In sharp contrast, the 
bandwidth-control Mott transition shows a clear first- 
order character with a jump of the double occupancy 
(for example at U c = 3.25 ± 0.05 for t' = -0.2). 

The phase boundary of metals and the Mott insula- 
tor shows that the charge gap A c opens at U = U c 
and shows marked linear dependence on U for U > U c , 
namely A c ~ a(U — U c ). The resultant V-shaped sin- 
gular phase boundary in the plane of U/t and [i (as 
seen in Fig. 11) is shown to be consistent with the sharp 
contrast between the continuous filling-control transition 
and the first-order character of the bandwidth-control 
transition. A general relation of the slope of the metal- 
insulator transition line in the [i-U phase diagram and 
physical quantities are also derived: In the case of the 
first-order metal- insulator transition, SU/S/j is expressed 
by the ratio of the jumps in the filling and the double 
occupancy An/ AD. In the case of the continuous tran- 
sition, 8U /S/j, is expressed by the ratio of the compress- 
ibility and dn/dU in the metallic phase. These relations 
support that the V-shaped phase boundary is resulted 
from the first-order BCMT coexisting with continuous 
FCMT with diverging compressibility. Experimental re- 
sults of Ti perovskite compounds are favorably compared 
with this V shape and the contrast between BCMT and 
FCMT. 
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Appendix A: BCS wavefunctions in the GPIRG 
framework 

In this appendix, we show the BCS wavefunction in 
the GPIRG framework and make some comments. 
The BCS wavefunction is defined by 

I^bcs) - I] ( Uk + v * c U c -kl) l°>> (AT) 

k 

where Uk and Vk satisfy u\ + \vk\ 2 = 1. 

By the canonical transformation of eq. (2), the BCS 
wavefunction is transformed as follows: 15 ^ 

iv'bcs) = n («* + vkcUk) n 4 io) 

k k' 

= n(«*4+«*4)|0> (A-2) 
k 

The BCS wavefunction is given as the 2V x N Slater 
matrix: 



'\jk 



u k e" 



ik R, 



for j 
for j 



1,...,V 
N+1.....2N, 



where k specifies N points in the Brillouine zone. Here 
Uk and Vk have the following forms: 



u k = 



where = Ek — H with Sk = — 2t[cos(k x ) + cos(fc y )] + 
At' cos(k x ) cos(ky) and [i is the chemical potential. The 
gap function can be given for several symmetries: For 
example, 




A fe = Ax < 



1, 



isotropic s wave 
s wave (cross) 
d x 2_ y 2 wave 
s wave (diagonal) 
d xy wave 



(cos(fcr) + COs(ky)), 

(cos(k x ) - cos(ky)), 
2 cos(k x ) cos(fcj / ), 
2 sin^) sin(fc y ). 

By taking the superposition of the columns of [4>]jk, [<j>]jk 
can be transformed into the real one: For example, in 
case that two columns are specified by k and k, they 
are expressed by cos(k • Rj) and sin(k • Rj) by using the 
crystal inversion symmetry. 

In the Hubbard model with on-site Coulomb repulsion, 
we note that there is no BCS-type mean-field solution. 
Hence, the most optimal state for L = 1 is not provided 
by the Slater determinant of eq. (AT). 

Appendix B: A proof of the Stratnovich- 
Hubbard transformation which 
hybridizes c and d particles 

In this appendix, a proof of the Stratnovich-Hubbard 
transformation which hybridizes c and d particles, eq. (9) 
is given. 

Let us rewrite eq. (9) with the site index omitted: 



exp 



-A T U | i (c^ + dU) - cWdj 
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= \ E exp[i/3 S ( C trf + dt c )] , (B-l) 



«=±i 



with (3 = cos- 1 (exp(-A r C//2)) for A T U > 0. 

The left-hand site of cq. (B-l) is expanded by noting 
that c^c, Sd, and c^cSd commute each other, as fol- 
lows: 



exp 
exp 



-A T U 1 1 (c+c + dU) - cWd j 



A T E/ 



c^c I exp 



A T C/ 



dtd) exp (A T [/cWd) 



{l + (c- A ^/ 2 - l) c t c } {l + (c- A ^/ 2 - l)rft d } 



{i + 



l)cWi} 



1 + ( C -^ u ' 2 - 1) ( c t c + Sd - 2cWd) . 



(B-2) 



The right hand site of eq. (B-l) is expanded by noting the 
fact that (cU + d^c) 2n = + dU - 2c Wd and (cU + 
d)c) 2n+1 = c^d + d)c for n = 1, 2, oo, as follows: 

exp [i/3s (ctd + d+c)] 

=t+E^(^+^r> 



n=l 



= i + i £(-iy 



(J3s 



>2n+l 



n=0 



(2n+ 1)! 



( C td + dt c ) 7 



+E(-d 



„(^) 2 



(cWc^d- 2 C t c dtd) , 



(2n)! 

= 1 + isin(/3s) (c f d + d f c) 

+ {cos(/3s) - 1} (^c + cftd - 2c*cdU) . 

Then, after the summation of the Stratnovich-Hubbard 
variables s = 1 and s = — 1, we obtain: 



i £ exp[i/3 S ( C t c + dtd)] 



s=±l 



= 1 + (cos (3-1) (^c + d f d - 2 C t C dtd) . (B-3) 

To make eq. (B-2) and eq. (B-3) equal, the relation 

A T lP 

cos (3 — 1 = exp [ T - — ) — 1 , 



should be satisfied for A T ?7 > 0. This determines the 
value of (3 as (3 — cos" 1 (exp(— A T U /2)). By substitut- 
ing this (3 to eq. (B-2) and eq. (B-3), we finally obtain 
eq (B-l). 
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